NASA 

Technical 

Paper 

i 2323 

i'.'O./ 

May 1984 



AUG 2 8 1984 

MAY 2 2 198? 


Finite-Difference Fluid Dynamics 
Computer Mathematical Models 
for the Design and Interpretation 
of Experiments for Space Flight 


Glyn O. Roberts, 
William W. FowKs, 
and TimothyL. Miller 



^ & Air Fom 
aedc ubiww •> 
E4060M1-C-0W4 







NASA 

Technical 

Paper 

2323 

I 1984 


Finite-Difference Fluid Dynamics 
Computer Mathematical Models 
for the Design and Interpretation 
of Experiments for Space Flight 


Glyn O. Roberts 

Roberts Associates, Inc. 

Vienna, Virginia 

William W. Fowlis 
and Timothy L. Miller 

George C. Marshall Space Flight Center 
Marshall Space Flight Center, Alabama 


NASA 

National Aeronautics 
and Space Administration 


Scientific and Technical 
Information Branch 


ACKNOWLEDGMENTS 


The authors wish to acknowledge support from the Universities Space Research Association at 
Columbia, Maryland, which was funded through contracts NAS8-34008 and NAS8-34315 from NASA 
Marshall Space Flight Center, Alabama. Also they acknowledge support from MSFC contract NAS8- 
34715 with Science Applications, Inc., of McLean, Virginia. This research was supported by the Global- 
Scale Atmospheric Processes Research Program of the Office of Space Science and Applications, and by 
the Microgravity and Science Applications Division, NASA Headquarters, Washington, D.C. 


iii 


TABLE OF CONTENTS 


Page 

I. INTRODUCTION 1 

II. THE ATMOSPHERIC GENERAL CIRCULATION EXPERIMENT 1 

A. Experiment 1 

B. Numerical Models 2 

C. Results 4 

III. CONVECTION IN A FLOAT ZONE 5 

A. The Float Zone System 5 

B. The Numerical Model 5 

C. Results 5 

IV. THE BRIDGMAN-STOCKBARGER CRYSTAL GROWING SYSTEM 7 

A. The Bridgman-Stockbarger Configuration 7 

B. Numerical Model 7 

REFERENCES 9 



LIST OF ILLUSTRATIONS 


Figure Title Page 

1 . Schematic drawing of the proposed AGCE apparatus 2 

2. A comparison of experimental and numerical results for the regime diagram for 

the baroclinic cylindrical annulus flows taken from Miller and Gall 3 

3. Preliminary regime diagram for the hemispherical-layer annulus configuration 4 

4. Flow in a rotating float zone driven only by surface tension variation 

in the outer cylindrical boundary 6 



TECHNICAL PAPER 


FINITE-DIFFERENCE FLUID DYNAMICS COMPUTER MATHEMATICAL MODELS FOR THE 
DESIGN AND INTERPRETATION OF EXPERIMENTS FOR SPACE FLIGHT 


I. INTRODUCTION 


Spacelab flights using NASA’s Space Shuttle have begun and regular flights are planned for the 
future. Among other uses, Spacelab has been designed to exploit the microgravity environment of an 
orbiting vehicle for science and technology. In particular, there are many basic fluid dynamics experi- 
ments, and materials processing studies involving fluid motions, which can only achieve their full potential 
in a low-gravity environment. NASA has also proposed to build a space station for which low-gravity 
experimentation will be a major application. This paper deals with three fluid dynamical experiments 
which are under consideration for space flight and with the computer mathematical models which are 
being developed and used for the design and interpretation of these experiments. 

Recent advances in numerical modeling methods and the continued improvement in computing 
speed and memory are opening up a new era for theoretical physics and engineering. In particular, 
accurate numerical models of non-turbulent, incompressible fluid flows in simple geometries can now be 
developed in a systematic manner. Sufficient numerical methods and numerical stability criteria are 
known so that this is no longer an art. Further, when accurate and detailed data are required for a fluid 
dynamical problem, numerical modeling is often the only procedure available. Nonseparability and non- 
linearity of the governing differential equations can make it impossible to use analytical methods, and 
the properties of many fluids and thin boundary layers can make it very difficult to use experimental 
methods. When such difficulties arise, the following procedure has been used to solve the problem [1,2]. 
First, develop a computer code and build an apparatus for the flow to be investigated. Then check the 
code predictions using data from relatively simple but significant measurements made with the apparatus. 
To obtain good agreement, it may be necessary to adjust some of the code control parameters, such as 
the space and time resolutions. Once validated, the code alone can be used to obtain accurate and 
detailed data about the flow over the range of validation. (There may also be occasions when a code can 
be validated using theoretical results over their range of validity.) 

The constraints and high cost of space flight experimentation mean that quantitative design 
studies should be performed. The procedure outlined previously is well suited to the design and inter- 
pretation of such experiments. Clearly, the validation method cannot now be carried out directly. 
However, results from one or more similar laboratory experiments which deal with the physical processes 
involved can be used instead to check the code. This procedure is being used to design a spherical baro- 
clinic flow model experiment of the large-scale atmosphere flow for Spacelab and to study the processes 
at work in crystal growing systems which are also candidates for space flight. 


II. THE ATMOSPHERIC GENERAL CIRCULATION EXPERIMENT 

A. Experiment 

A substantial amount of research into large-scale baroclinic atmospheric dynamics has been per- 
formed using a laboratory cylindrical annulus model [3], The microgravity environment of Spacelab 
presents the opportunity to realize a true spherical-layer model. Such an experiment, known as the 



Atmospheric General Circulation Experiment (AGCE), has been proposed to NASA for Spacelab flights. 
In the AGCE radial gravity is simulated by a dielectric body force, which is only dominant in the absence 
of ordinary gravity. The AGCE has been described elsewhere [4,5,6] and a substantial number of scien- 
tific support studies have been performed [7-13], A spherical convection experiment known as the 
Geophysical Fluid Flow Cell which also uses the radial dielectric body force will be flown on Spacelab 3 
[14], Proposed configurations for the AGCE are shown in Figure 1. For quantitative design studies 
numerical models are required, and a numerical design studies program is well underway [15,6]. The 
most important design criterion is that baroclinic instability be realized in the apparatus. 



Figure 1. Schematic drawing of the proposed AGCE apparatus. The high latitude boundary 
will be used if a hemispherical annulus configuration is selected. 


B. Numerical Models 

To determine in parameter space the baroclinically unstable region for the AGCE, two finite- 
difference codes have been developed. The first solves the axisymmetric Navier-Stokes equations in 
spherical coordinates to obtain basic states. The second uses the linearized, three-dimensional Navier- 
Stokes equations to examine these basic states for instabilities with azimuthal structure. Detailed descrip- 
tions of these two codes have been given by Fowlis and Roberts [6], Only a very brief description is given 
here. 


With appropriate choices for the geometrical control parameters, we can calculate basic state 
flows in a spherical or hemispherical layer, or in a sphere. Boundaries can be introduced at any latitude. 
We can also approximate a cylinder or a cylindrical annulus by choosing a very large inner radius and a 
limited latitudinal range from the axis of symmetry for the latitudinal boundaries. A wide range of 
boundary conditions (e.g., no flow, stress free, insulating, and conducting) is available. Nonuniform 
meshes and implicit iterative methods are used. The iterative method for steady solutions is based on 
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time-stepping but has the options of different time steps for velocity and temperature and of a time step 
varying smoothly with position according to specified powers of the mesh spacings. This allows for more 
rapid convergence. The linear azimuthal wave stability equations are solved by an iterative method which 
corresponds closely with the axisymmetric method. 

Validation of the axisymmetric code using data from previous work in cylindrical geometry was 
described in Reference 6. To check the spherical terms we compared the code predictions for spin-up 
in a homogeneous sphere of fluid [16] with laser-Doppler measurements of this flow [17]. Good agree- 
ment was obtained. To validate the stability code, it was applied to the computations of basic state flows 
in a rotating cylindrical annulus. Good agreement was obtained with the experimental results of Fowlis 
and Hide [18]. A parallel research effort, in which a similar code was used and approximately the same 
agreement obtained, has been published by Miller and Gall [13]; their results are shown in Figure 2. 



Figure 2. A comparison of experimental and numerical results for the regime diagram for the 
baroclinic cylindrical annulus flows taken from Miller and Gall [13]. The experimental 
results for marginal stability are shown by the continuous curve and the numerical 
results by the dashed curve. Further details can be found in the paper. 
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In further support of the AGCE program, a fully nonlinear, three-dimensional code is being 
developed. This code uses the same spherical-layer geometry and allows for the same geometrical options 
as before. It is designed to compute either steady states, using the same techniques as before for rapid 
convergence, or time-dependent solutions. The primitive variables are used on a nonuniform, staggered 
computational mesh. Advection, diffusion, the Coriolis force, and internal waves are all handled impli- 
citly, to reduce or eliminate the corresponding limitations on the time step. This code will be validated 
using known wave flows in the cylindrical annulus. It is important to have prior estimates of the ampli- 
tude and structure of the wave flows to be expected in the AGCE. The code will be used both for 
design of the AGCE and later for interpreting actual AGCE flows. 


C. Results 

A definite configuration for the AGCE in terms of latitudinal boundaries, radial depth, boundary 
temperature distributions, and fluid properties is not yet settled; such specifications should come out of 
the design studies. Figure 1 shows a hemispherical-layer configuration, which can be heated and cooled 
at the equatorial and high latitude boundaries, respectively, similar to the cylindrical annulus. We are 
computing the regime diagram for such a spherical annulus configuration, and the results to date are 
shown in Figure 3. At this point, the results suggest that the diagram will be very different from the 
cylindrical annulus diagram shown in Figure 2. 



Figure 3. Preliminary regime diagram for the hemispherical-layer annulus configuration. The azimuthal 
wavenumbers are shown, and whether they were found to be stable (S) or unstable (U) is also 

indicated. The results are shown as a plot of agDAT/H“L^) versus ir^(= 

for a fixed value of Pr (= v/k), where D is the radial fluid depth, L the latitudinal extent, a 
the coefficient of cubical expansion, v the kinematic viscosity, k the thermal diffusivity, 

AT the imposed latitudinal temperature difference, and £2 the rotation rate. For these 
results the inner radius = 4.5 cm, the outer radius = 6.5 cm, and the high latitude 
boundary was located at 0.45 radians from the polar axis. The liquid was 

silicone fluid with v = 0.01 cm^/s. 
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Problems are being experienced with slow convergence of the algorithm which computes the 
fastest growing (or slowest decaying) eigenmode for linear disturbances to the basic state. The algorithm 
is related to time-stepping the evolution equations, but with a time step which varies with position and 
for different variables. This attempt to obtain rapid convergence has been successful for a large number 
of cases, particularly in an annular geometry. However, it has not been as successful for most of the 
spherical-layer geometry cases analyzed. Therefore, small time steps with correspondingly slow 
convergence must be used. 


III. CONVECTION IN A FLOAT ZONE 

A. The Float Zone System 

In the float zone process of crystal growth, a ring heater maintains a small zone of liquid melt 
between a polycrystalline feed rod and a rod of single crystal. The small volume of liquid is held in 
place by surface tension. The feed rod is moved into the zone where it melts, while the other rod is 
pulled out of the zone resulting in growth as a single crystal. In order to achieve desired properties for 
the crystal, the liquid melt is usually made up of two or more different substances. Ideally, the com- ^ 

ponents should be distributed uniformly throughout the crystal. However, in general, this is very diffi- 
cult to achieve [19]. The variation of component concentration in the crystal is caused by the flow and 
temperature fields in the melt. Strong flows can also cause crystal imperfections with a single component. 

Since buoyancy-driven flows are suppressed in the absence of gravity, float zone systems are candidates 
for space flight. However, even in the absence of gravity, the variation of surface tension with tempera- 
ture on the free surface will drive motions. 

As a first step in understanding the interaction between the flow and temperature fields and the 
crystallization process, Smith and Greenspan [20] investigated the motions occurring in the float zone 
in the absence of crystal growth. These workers added rotation to examine its influence in confining 
the flows to boundary layers. Adding rotation also adds another driving force, centrifugal buoyancy. 

Smith and Greenspan examined, separately, the flows driven in a rotating cylinder of fluid by thermo- 
capillary and centrifugal buoyancy forces. Ordinary gravity was neglected. These workers used analytical 
methods and linearized equations. They considered a melt zone which was hot at one end and cold at 
the other. Using numerical methods, this model is extended to the nonlinear case; and, using both 
analytical and numerical methods, the more realistic model with a temperature distribution symmetrical 
about mid-depth is examined. 


B. The Numerical Model 

For this work, the axisymmetric code already developed for the AGCE basic state studies can be 
used (see Section II). The conversion to cylindrical geometry was again carried out by letting the radius 
of the inner sphere be very large and considering only a small latitude range from the axis of symmetry. 
The stress boundary conditions at the outer cylindrical surface were changed to implement the thermo- 
capillary boundary condition (temperature-dependent surface tension). 


C. Results 

Following Smith and Greenspan [20], thermocapillary-driven flow and the centrifugal buoyancy- 
driven flow were considered separately. Smith and Greenspan presented results for a cylinder of silicon 
liquid (for details see Figure 4 caption). They found, for the linearized problem, that the flow driven 
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by the temperature variation of surface tension at the free outer cylindrical surface is confined to a 
boundary layer on that surface and that there is no flow in the interior. However, Smith and Greenspan 
noted that the dimensionless parameter determining the strength of the nonlinear terms was not small 
for the case they considered. Their case was repeated using the fully nonlinear numerical model and it 
was found that the boundary layer flow forces its way into the interior. The stream function is shown 
in Figure 4. Additional computation showed that for this case there is little difference between the 
rotating and the nonrotating flows. However, the suggestion to use rotation to limit the extent of the 
flow may still be useful if a suitable surfactant, to reduce the thermocapillary effect, can be found. 

This paper does not present the rest of the results on float zone convection. These results should 
be published soon [21]. 

AXIS 

OF 

SYMMETRY 



Figure 4. Flow in a rotating float zone driven only by surface tension variation in the outer 
cylindrical boundary. Contours of the stream function are plotted. The results are for: 
radius = 1 cm, height = 1 cm, vertical temperature difference = 5°C, and rotation = 

1 rad/s. The fluid was liquid silicon with surface tension = 720 dyne/cm, surface 
tension temperature coefficient = -0.43 dyne/cm C°, kinematic viscosity = 

0.0035 cm“/s, and thermal diffusivity = 0.15 cm^/s. The maximum value 
for the vertical (axial) flow at the surface was about 25 cm/s. 
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IV. THE BRIDGMAN-STOCKBARGER CRYSTAL GROWING SYSTEM 


A. The Bridgman-Stockbarger Configuration 

The Bridgman-Stockbarger system is one of the most widely used methods for growing crystals 
in the laboratory. In this method, a cylindrical ampoule containing the liquid melt is lowered steadily 
from a hot isothermal container into a relatively cooler container. Solidification begins at the bottom 
of the ampoule and proceeds upwards as the ampoule is withdrawn. Again as for the float zone system, 
to achieve desired properties for the crystal, the liquid usually consists of two or more different sub- 
stances. Although the Bridgman-Stockbarger system has been used successfully to produce crystals for 
many applications, there are also many cases where the full potential of the crystalline material is not 
achieved for the reasons already given in Section III. 

The temperature and flow fields influence the crystal growth in two ways. First, they can affect 
the shape of the solidification interface which in turn can affect the radial concentration of the melt 
components. Second, the flow field can modify the thickness of the diffusion boundary layers of com- 
ponent-rich liquid which form near the solidification interface, as a result of differential segregation of 
the components during solidification. This thickness modification then affects the concentrations of the 
melt components incorporated into the crystal. 

A fuller understanding of the physical processes at work in this system would almost certainly 
lead to improvements in the crystals produced, and a numerical modeling effort appears to be a worth- 
while approach. For many cases of practical interest, experimental measurement and observation are 
difficult or impossible, since the substances are metals at elevated temperatures and pressures. Analytical 
approaches are hindered by the complex geometry and by the nonlinearity and instability of the flows. 


B. Numerical Model 

A code is being developed to compute time-dependent and steady state axisymmetric solutions 
for the Bridgman-Stockbarger system. The code will be flexible, to describe a variety of geometries, 
materials, and thermal configurations, for laboratory and industrial applications and for space. 

The temperature equation is solved in the liquid, the solid, and the ampoule, using flexibly- 
defined thermal forcing on the boundaries and using appropriate interface conditions including latent 
heat release. The concentration equation for a secondary component and the equations of motion are 
solved in the melt. An incompressible Boussinesq fluid, with the buoyancy in general a function of both 
the temperature and concentration, and with an imposed density (and volume) change on solidification 
is assumed. 

The code uses a nonuniform mesh defined using a coordinate transformation, so that the 
unknown interface is a coordinate surface. The radial mesh is also nonuniform. The interface height is a 
function of radius and time; its evolution is determined using an implicit representation of the interface 
conditions, based on the equations of state. 

The heat equation is time-stepped using an implicit method. For typical applications, the thermal 
conductivities are very large, and this implicit method essentially amounts to a separate solution of the 
steady state problem at every time step. The concentration equation is also time-stepped by an implicit 
method. Here the diffusivity is very low, and this method is normally required only if the flow is rapid 
and the advection stability criterion would otherwise be violated. 
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The equations of motion are solved using a vorticity and stream function formulation, with 
viscosity, advection, Coriolis force, and internal waves all treated implicitly. Again, this is required 
because the time step would otherwise have to be very small. The Bridgman-Stockbarger system will 
require prolonged calculations. Model predictions will be compared with laboratory measurements made 
using different Bridgman-Stockbarger systems by our colleagues in the Space Science Laboratory at 
Marshall Space Flight Center (S. L. Lehoczky, D. O. Frazier, F. R. Szofran, and D. Chandra). 
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